Feedback-induced oscillations in one-dimensional colloidal transport 
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We investigate a driven, one-dimensional system of colloidal particles in a periodically currogated 
narrow channel subject to a time-delayed feedback control. Our goal is to identify conditions 
under which the control induces oscillatory, time-periodic states. The investigations are based on 
the Fokker-Planck equation involving the density distribution of the system. First, by using the 
numerical continuation technique, we determine the linear stability of a stationary density. Second, 
the nonlinear regimes are analyzed by studying numerically the temporal evolution of the first 
moment of the density distribution. In this way we construct a bifurcation diagram revealing the 
nature of the instability. Apart from the case of a system with periodic boundary conditions, we 
also consider a microchannel of finite length. Finally, we study the influence of (repulsive) particle 
interactions based on Dynamical Density Functional Theory (DDFT). 



I. INTRODUCTION 

The study of transport of particles in complex geome- 
tries is a major topic in nonequilibrium statistical physics 
with relevance in diverse fields such as biology, condensed 
matter and nanotechnology P, Q ■ Exemplary systems 
are colloids in optical (or otherwise modulated) poten- 
tials (bio-)molecules in microchannels pj, cold 
atoms in optical lattices Q, and magnetic particles ad- 
sorbed on ferrimagnetic garnet films Q. Depending on 
the details of the (often one-dimensional) potential, a va- 
riety of fascinating effects can been observed, including 
ratchet mechanisms (Toj , giant diffusion [TT| , and anoma- 
lous (subdiffusive) transport pT3 - fT5l ] - For colloids, which 
are typically of the size of nano- to micrometer, many of 
these effects can be monitored by real-space experiments 
(see, e.g., EMI). 

In the present paper we investigate a one-dimensional 
colloidal system where the force exerted by the (static) 
modulated potential is supplemented by a feedback con- 
trol force, i.e., a force depending on the state of the sys- 
tem. Feedback control in the context of Brownian sys- 
tems is a focus of growing interest, and the method has 
already been a ppl ied, on a theoretical level, to Brown- 
ian motors (20l. l2l| and flashing or rocking ratchets [22| . 
Moreover, a first experimental realization of a feedback- 
controlled flashing ratchet already exists [23| . The overall 
goal of the feedback control in this context is to manip- 
ulate and/or optimize transport properties such as the 
current in a flashing ratchet. Beyond this more appli- 
cational motivations, feedback-controlled transport phe- 
nomena are also of fundamental interest due to the sub- 
tle interplay between state-dependent control protocols, 
thermodynamics and information theory (24l . [25j . Very 
recently, a colloidal system under feedback control was 
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used as experimental realization of an "information heat 
engine" [26l | . 

Following an earlier study of two of us [13] , we here em- 
ploy a feedback control with delay where the control term 
involves the difference between an system variable (the 
control target) at time t and its value at time t — T, with 
r being the delay time. Such a time delay often occurs 
in experiments due to the time lag between measurement 
and feedback. In the general context of control of non- 
linear system [28|, time-dela yed feedback control (which 
was introduced by Pyragas [29| ) has been proven to be 
extremely efficient e.g. for the stabilization of chaotic or- 
bits; however, steady states can be manipulated as well. 
In [23] we used the delayed feedback control to manip- 
ulate the current of interacting colloids driven through 
in a tilted washboard potential. Indeed, it turned out 
that the control can optimize the current and even yield 
a current reversal, similar to what has been previously 
seen for non-interacting systems [30l l3l| . 

In the present study we go one step further and ask 
to which extent the time-delayed feedback control can 
induce dynamic states not seen in the uncontrolled sys- 
tem, which involves a purely static potential. Specifically, 
we search for the existence of spatio-temporal structures 
characterized by an oscillatory distribution of particles. 
That time-delayed feedback control itself can indeed gen- 
erate novel dynamics has recently also been seen in other 
extended systems [32|. 

Our investigations are based on the numerical solution 
of the nonlinear Fokker-Planck equation combined with 
a linear stability analysis. Thus, the basic dynamic vari- 
able is the time-dependent density field. Similar to our 
previous study p7| our control protocol involves the aver- 
age particle position as a control target, a quantity which 
is, in principle, accessible by experiments. We investi- 
gate both, infinite systems, that is, colloidal particles on 
a ring (see [H| for an experiment), and a finite system 
to which we refer as microchannel geometry (33|. For 
both geometries, we do indeed find oscillatory states at 
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appropriate system parameters and finite values of the 
control strength and the delay time. The character of 
these oscillations, on the other hand, strongly depends 
on the set-up. 

In the last part of the paper, we briefly discuss the im- 
pact of repulsive interactions between the particles. This 
is done on the basis of Dynamical Density Functional 
Theory (DDFT), where the microscopic interactions en- 
ter via the free energy functional. Indeed, in the last 
years DDFT has been applied to a variety of driven col- 
loidal systems J34[ , including attracting colloidal particles 
in ID (time-dependent) ratchet potentials [35j . 

The rest of the paper is organized as follows. Sec. [TT1 
contains the formulation of the problem. In Sec. IIII| 
we determine the stability of a nontrivial stationary, 
spatially-periodic density distribution by linearizing the 
Fokkcr-Planck equation subject to a time-delayed control 
force. For the nonlinear regimes, we provide numerical so- 
lutions of the (full) Fokker-Planck equation and construct 
a bifurcation diagram based on monitoring the amplitude 
of the solutions. In Sec. HVl we discuss the impact of the 
channel geometry on the oscillation instability. For that 
purpose, we consider a microchanncl of finite length L 
and follow the same procedure as for the infinite system. 
The role of interactions between the particles is briefly 
discussed in Sec. |V] 

II. THEORY 

Our model system consists of overdamped colloidal par- 
ticles in a one-dimensional channel of length L. The parti- 
cles are subject to a spatially periodic, symmetric "wash- 
board" potential i/ w b(z) = Uq cos 2 (kz), where k defines 
the wavelength and Uq is the amplitude. Henceforth we 
set ka = 1 where a corresponds to the effective (gyration) 
radius of the particle. The washboard potential is tilted 
by a constant force Fbias = Fqz (with z being the unit 
vector in z-direction) corresponding to a linear potential 
Ubms = —Fqz. The tilting leads to an effective motion of 
the particles in the direction of sign(i 7 o) along the z-axis. 
In addition to these static potentials, we assume that the 
particles are subject to a time-delayed feedback control 
force of the form [27| 

F lh (t, t) = -K (l - tanh [f(t) - f(t - r)] ). (1) 

where f(t) is a space- averaged coupling function which 
depends on the internal dynamical variables of the sys- 
tem. Specifically, we set 

/(*)= / f(z)p(z,t)dz, (2) 

J Zq 

with zq being the coordinate of the origin of the z-axis. 
Note that when the feedback term is switched on, the 
effective constant driving force is given by 7 = F — Kq. 

In Ref. lU we used a force of the type ([T]) to manip- 
ulate, or, more precisely, to revert the net current in- 
duced by the biasing force. Here, we aim to explore to 



which extent such a feedback force can induce oscillatory 
(time-periodic) density states, which correspond to syn- 
chronized oscillations of the particles along the channel. 
We note that the control force in Eq. ([TJ is of Pyragas 
type |2{|, i.e. the difference of the "control target" at 
time t and its value at time t — r is used as input for 
the feedback loop, where r is called "delay time" and Kq 
is the control amplitude. We stress that the ansatz for 
the feedback force in Eq. ([1) is clearly of heuristic nature, 
and thus cannot be derived from any physical potential. 

Collecting all external contributions gives the total ex- 
ternal potential 

U cx t =U wh (z) + Ufb(t, T) + U hiaS (z) 

=U cos 2 (z) + K z(l - tanh [f(t) - /(* - r)]) - F z, 

(3) 

where we assumed that the contribution of the feedback 
force is linear in the position coordinate z. In the absence 
of interactions between the particles, their Brownian mo- 
tion is governed by the following Fokkcr-Planck equation 

m 

p-i dp(z, t) d 2 p(z,t) d 

r Qt =k B T—^- 2 — \p{z,t)ft{z,t;T)], (4) 

where the drift coefficient ft can be calculated from the 
external potential f7 0Xt via ft = —U' cxt where ' denotes the 
derivative w.r.t. z. The mobility coefficient in Eq. (U) is 
related to the diffusion constant via V = /3Dq [where (3 = 
I/(fceT)] and we set its value to Fre^T/cr 2 = 1. Time 
is measured in units of the Brownian time scale tb = 
cr 2 /(rfcsT), which is of the order of I0 _9 s for typical 
Brownian particles. 

In the following, we choose a positive value for Fq such 
that the particles move preferentially to the right when 
the control force is being switched off. Also, for the main 
part of our investigations, we consider the system to be 
infinitely extended ("circular ring"), that is, we apply 
periodic boundary conditions over a length L, which is 
a multiple of a period of the washboard potential, i.e. 
L = mra, with n = 1,2,.... However, in Sec. IIV1 we 
also discuss the impact of finite channel length. 

III. DELAY-INDUCED INSTABILITY 

In the absence of control (Kq = 0), the driven system 
(F > 0) settles into a stationary, non-oscillatory state. 
The corresponding distribution p s (z) can be found ana- 
lytically [K|. In the presence of control, the behavior 
of the system depends on the interplay of the parame- 
ters of the control term, on the one hand, and the tilted 
washboard potential, on the other hand. For certain pa- 
rameter combinations we still find a stationary state, in 
which the difference f(t) — f(t — r) disappears. How- 
ever, as we will demonstrate, for suitable parameters the 
stationary distribution may become unstable, leading to 
a stable time-periodic distribution p(z,t). To interpret 
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this instability, it is useful to reconsider Eqs. (Q][2]) from 
a somewhat different perspective. In particular, since the 
control target fit) involves a spatial integral over the en- 
tire distribution p(z,t), we have effectively introduced a 
coupling between the colloidal particles. Moreover, this 
coupling is of infinite range (within the periodic system 
considered). Thus, the time-delayed feedback control in- 
troduces an effective interaction of mean-field type, and 
it is this interaction, which may lead to new stationary 
states as well as new dynamic regimes, associated with 
time-periodic density oscillations. We note that the par- 
ticular type of feedback used here is of purely heuristic 
nature. This implies that the Fokker-Planck equation 
Eq. (TJJ cannot be directly translated into a correspond- 
ing system of coupled Langevin equations. 



A. Linear stability analysis 

We start by applying a linear stability analysis in or- 
der to investigate the impact of the delayed feedback con- 
trol. To this end, we rewrite the dimensionlcss Fokker- 
Planck equation in terms of the effective potential U e s = 



Uq cos z + (Kq — Fq)z as follows 

dp{z,t) d 2 p(z,t) d 
dz 



ot 



dz 2 



p(z,t) 



dU cff 



Xotanh [/(t)_/(t_r)] 



dp(z,t) 
dz 



The non-trivial stationary state p s (z) satisfies 







d_ 

dz 



da 



err 



dp s (z) 



dz dz 

which can be written as an eigenvalue problem 
Lp s = Xp s , 



(5) 



(6) 



(7) 



with the stationary Fokker-Planck operator L = d 2 + 
(d z U e e)d z + d 2 U e g and zero eigenvalue A = 0. 

We are interested in the onset of an oscillatory insta- 
bility of p s . Following the approach, developed earlier 
[37l HH , we set 

p(z,t) = p s (z) + e(C(z) cos ujt + S(z) smut) , (8) 

where e is the (small) amplitude of the perturbation, u is 
the unknown onset frequency and the unknown functions 
C{z) and S(z) determine the shape of the perturbation. 
With this ansatz, the Eq. (|5|) can be linearized in e to 
yield 

— ljC sin wt + ujS cos tut =C" cos uit + S" sin out 

+d z [U' eS (C cos ut + Ssinut)) 
-K oP ' s [fit) - /(* - r)] , (9) 

where / stands for the derivative w.r.t. z and the per- 
turbed mean field is given by 



fit) = (C) cosujt + (S) smut, 



(10) 



with (C) 

rL+z 



= J^ Za f(z)C(z)dz and (S) = 

Iz +Z ° f( z )S( z ) dz. The initial moment of time can 
always be chosen in such a way that, for instance, 
fit) = cosujt. This implies two additional integral 
conditions on the functions C(z) and S(z) 



L+z 



f(z)C(z)dz = l, 



L+zq 



f(z)S(z)dz = 0. (11) 



Then the difference f(t) — f(t — r) becomes 

f(t) — fit — r) = cosoj<(1 — coswt) — sin ut sin lot. 

(12) 

Finally, equating the coefficients of sinwi and coswi in 
Eq. ([9]) , we obtain two coupled equations for the unknown 
functions C[z) and S(z) 

S" =-u>C-d z [U cfi S} - K smuT P ' s 

C" =luS - d z [U cS C] + K {1 - cosujt)p' s . (13) 

In order to find the stability threshold, one needs to 
solve Eqs. (fT3"|) simultaneously with Eq. (J7J. To this end 
we proceed as follows. We rewrite Eqs. (|7I13|) as an au- 
tonomous dynamical system of seven first order equa- 
tions, including the equation for z, which takes the form 
z' = 1. The total number of the system parameters is 
thereby extended by two additional parameters, namely 
by the onset frequency u> and the (zero) eigenvalue A. 
The above dynamical system of seven equations is sup- 
plemented with three integral conditions on the functions 
p s (z), C(z) and S(z). These are given by Eqs. (jTTj) and 
by the normalization condition on p s 



L+z 



p s (z) dz = N, 



(14) 



where N is the normalization parameter. The boundary 
conditions for all involved functions are taken to be peri- 
odic in the interval z € [zq, zq + L\. This boundary value 
problem (BVP) of seven equations and three integral con- 
ditions is then solved using the numerical continuation 
technique (AUTO) [11] (see AppendixlAl for details). 



B. Stability thresholds 

Before proceeding, it is important to notice that the 
linear stability of the stationary distribution p s crucially 
depends on the choice of the coupling function f(z) [see 
Eq.©]. Thus, the solution of the BVP Eqs. (JT3J) is invari- 
ant under the shift of the coordinate system z — > z + 5, 
with arbitrary 5, only if the coupling function f(z) is 
itself periodic with the period L. Indeed, the inte- 



grals f*° +L f(z)C(z) dz and f*° +n f(z)S(z) dz, are shift- 
invariant if f{z) is L-periodic. On the contrary, if the pe- 
riod of f(z) is different from L, or if f(z) is a non-periodic 
function, then the stability threshold depends on the par- 
ticular choice of origin of the z-axis, i.e. it depends on 
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Figure 1. (Color online) (a) Stability threshold in the plane 
(7, Uo), obtained for r = 1 and different values of Ko, as given 
in the legend. The stationary state is unstable in the area en- 
closed by the respective curves (shaded area for Ko = 7). (b) 
Stability threshold in the plane (Ko, Uo) for zero drive 7 = 
and different delay times as in the legend. The instability re- 
gion always lies to the right from the respective curve (shaded 
area for r = 0.1) 



zq. Following |27j we use a linear, non-periodic coupling 
function f(z) = z. For further calculations, the origin of 
the z-axis is chosen in the maximum of the washboard 
potential U(z), implying that zq = 0. 

First, we fix t = 1 and compute the stability threshold 
in the plane of parameters (7, Uo), where 7 = Fq—Kq = 0, 
for three different values of the coupling strength Ko = 
7, 10, 12, as shown in Fig.QJa). The stationary density, 
normalized with N = 1, i.e. J Q p s (z)dz = N = 1, is 
unstable in the regions bounded by the corresponding 
closed curves. Thus, for K = 7, the instability occurs in 
the shaded area. As expected, the area of the instability 
expands if the coupling strength K is increased. 

Interestingly, the stability diagram for the case when 
the origin of the z-axis is chosen in the minimum of the 
washboard potential (zq = 7r/2), can be obtained from 
Fig. (Ha) by the transformation Uo —> —Uo- For any other 
choice of z$, the topology of the stability threshold is 
much more complex and generally contains four different 
bounded regions (not shown). 

The effect of the time delay is demonstrated in Fig.[IJb) 
for the choice 7 = 0. The stationary density is unstable 
in the area, which always stretches towards larger values 
of the coupling strength Ko, as shown by the shaded 
area for r = 0.1. Decreasing r leads to the suppression 
of the instability, which clearly demonstrates that the 
instability is induced by the presence of the time delay 
in the coupling term. 

However, it should be emphasized that having only a 
time-delayed coupling does not suffies to induce the in- 
stability. Thus, from Fig.[IJb) it follows that even if r 
and K are rather large, e.g. r = 200 and K ~ 100, 
the stationary density is linearly stable for a vanishingly 
weak or an infinitely strong washboard potential Uq- At 
fixed Kg, only a certain combination of r and Uo ren- 



Figure 2. (Color online) Three-dimensional view of Fig.QJa). 
The onset frequency ui as a function of Uo and 7 for r = 1 
and three different Ko- The projection onto the (7, Uo) plane 
recovers Fig.[]Ja). 



ders the system unstable with respect to an oscillatory 
perturbation. Consequently, the onset of the synchro- 
nized time-periodic state is a the effect of the combined 
action of the time-delayed coupling and stationary peri- 
odic external modulation in the form of the washboard 
potential. 

By following the stability threshold in the parameter 
space, we additionally obtain the onset frequency uj di- 
rectly on the threshold. The latter carries an important 
information about the time scale of the newly born oscil- 
latory states, given by T — 2ir /ui. The three-dimensional 
view of Fig. [TJa), extended by the onset frequency uj, is 
shown in Fig. [5J It can be seen that the smallest tempo- 
ral period T ~ 1 corresponds to positive Uq, whereas the 
period of the time-periodic states born at negative Uo is 
up to three fold larger, i.e. T ~ 6. 



C. Nonlinear regime: Numerical study of the 
one-body distribution 

The observation of the linear instability with respect 
to oscillatory perturbations predicts possible deviations 
from the stationary state for a given ( "overcritical" ) pa- 
rameter set. In order to study the full nonlinear dynam- 
ics, however, the full Fokker-Planck equation of the sys- 
tem has to be solved. To explore these nonlinear effect 
we solve Eq. (0]) numerically. Specifically, we employ a 
standard "Forward-Time Centered-Space" (FTCS) finite 
difference method [4(| and integrate Eq. Q starting from 
an inital distribution p(z,t = 0). 

As initial density distribution we choose the equilib- 
rium density distribution in a one-dimensional wash- 
board potential with periodic boundary conditions and 
the external bias, as well as the control, being switched 
off. This implies 



p(z,t = 0) = po exp 



-fiUa cos' 



(;) 



(15) 



5 





0.18 




0.16 




0.14 




0.12 


N 


0.1 


a. 


0.08 


e 


0.06 




0.04 




0.02 








A 

N 

V 




120 



Figure 3. (Color online) Results for a controlled system in the periodic regime, (a) Density distribution p(z,t) as a function 
of position for selected times to/TB =0 (red curve), ti/m = 16.9 (green-dashed curve) and £2/75 — 68.6 (blue-dashed curve), 
(b) Average particle position {z)t as a function of time. The parameters are Uo = 7fcsT, Fq — IksT/a, Ko = TknT/a, r — tb 
and L — Wna. 



where po ensures the normalization condition J Q p(z, t = 
0)dz = N = 1. We consider a fixed control amplitude 
A'o = lk^T/a and focus on parameter values near the 
stability threshold [see Fig. 1(a)]. Specifically, we con- 
sider the "balanced case" 7 = F — K = and wash- 
board amplitude values of Uq = 8ksT (linearly stable) 
and Uo = 7ksT (linearly unstable), respectively. Begin- 
ning with the latter case, we plot in Fig.^a) snapshots of 
the density distribution p(z, t) for three subsequent times. 
The initial distribution p(z,t = 0) is periodic in the posi- 
tion coordinate z with a (spatial) period that is equal to 
the valley-to-valley distance of the washboard potential, 
that is A w b = iter. As expected, the values for the den- 
sity distribution are increased at the position coordinates 
^valley _ ( wn ere i = [0, 1, 2, 3, ... ]) corresponding to 
the valley positions of the washboard potential. The spe- 
cific times ti, ti are chosen such that the appertaining 
density distributions (shown as the green-dashed curve 
and the blue-dashed curve in Fig.[3Ja), respectively) have 
a maximum displacement from the initital configuration. 
Inspecting the curves, we find that after a response time 
of roughly 5tb, the system settles indeed into a stable 
time-periodic density state, i.e. p(z,t + T) = p{z,t), 
where the distribution oscillates around a washboard min- 
imum position with a maximum displacement of approxi- 
mately 0.33cr. The appearance of such stable oscillations 
is consistent with the stability diagram in Fig. [TJa). We 
supplement the discussion by plotting in Fig. [3jb) the 
particle position averaged over one period of the poten- 
tial, that is, {z) t = J Awb p(z,t)dz. Clearly, (z) t oscillates 
as a function of time with the same frequency to as the 
frequency of the density oscillations. From now on, we 
therefore use the function (z) t to obtain the cycle time 



T = 2tt/uj. 

We now turn to the case Uo = 8/csT. For this parame- 
ter, the perturbating forces, specifically the constant tilt- 
ing force F blas and the feedback force F f b -, do not lead 
to an oscillating state. This is illustrated in Fig. @] It is 
seen that the oscillations at early times are being damped 
resulting in a stationary, non-periodic density for times 
t > 400tb. In Fig. SJa) we show the results for the sta- 
tionary density profile as a blue-dashed curve, as well 
as for two additional profiles, where the red and green- 
dashed curves represent the initial and one intermediate 
state, respectively. As can be seen from Fig.^b), the av- 
erage particle position (z) t , converges to a constant value 
(z)stat/c ~ 0. In other words, the density displacement 
related to the equilibrium position at each potential val- 
ley vanishes and, thus, the profiles for times t^jTB = 400 
and to/rB = coincide. 

We note that for both values of Uo considered (cf. 
Figs. OH]), the dynamics at early times is transient. In 
the oscillatory case, the transient regime lasts for about 
t rcs ~ 5tb. For the linear stable case, the transient re- 
sponse time can be much larger. For example, for the pa- 
rameter values that we considered in Fig.|3]the stationary 
state is not reached for times less than i lcs w 400tb. 

We have repeated the numerical calculations described 
above for a range of parameters Uo and the choice 7 = 
Fq — Ko = 0. In this way we can construct a bifurca- 
tion diagram characterizing the nature of the instabil- 
ity. As a measure of the instability, we use the oscil- 
lation amplitude of (z)t- Specifically, we obtain local 
extrema values ( z ) max / mm from the function (z) t and av- 
erage these over several periods. The results are summa- 
rized in Fig [5J At large values of Uq we only find one 
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Figure 4. (Color online) Results for a controlled system in the regime where the final state is stationary, (a) Density distribution 
p(z,t) as a function of position for selected times £o/tb = (red curve), ti/rs = 50.4 (green-dashed curve) and tijrB = 400.0 
(blue-dashed curve), (b) Average particle position (z)t as a function of time. The parameters are Uo = 8fcsT, Fq = IksT/a, 
Kq — 7ksT /a, t = tb and L = IQiro. 



stable attracting fixed point ( z ) max / mm = o. However, 
by decreasing the washboard amplitude Uo the station- 
ary state (characterized by ( z ) max / min = for all times 
t greater than the transient time) loses stability and sta- 
ble limit cycle oscillations occur (see, e.g., Fig. [3]). This 
happens in an essentially continuous manner, as Fig. [5] 
reveals. We thus conclude that, upon decreasing Uo be- 
low [/q w 7.95/cbT the system undergoes a supercritical 
Hopf bifurcation. For Uq < U§ all resulting trajectories 
perform limit cycle oscillations about the former station- 
ary state ( z ) max / mm = o. We note that all neighboring 
trajectories approach the limit cycle. Thus, for Uo < U§, 
the limit cycle is stable and the only attractor in the 
system. Furthermore, due to the spatially left-right sym- 
metry in the potential, the local extrema of {z) t appear 
in symmetrical pairs at ±|(z) max |. 
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Figure 5. Local extrema values of the function (z)t as a func- 
tion of the washboard amplitude Uo- The parameters are 
Fa = TkBT/a, Ko = 7kBT/a, r = tb, L = lOrnr. 



D. Cycle time for the oscillating density state 

An interesting question is to which extent the density 
oscillation frequency depends on the different system pa- 
rameters such as washboard potential amplitude Uo and 
time delay r. In Fig. [2] we have already shown results for 
the onset frequency based on the linear stability analysis. 
Here we present corresponding data obtained in the non- 
linear regime. To this end we define the cycle time T as 
the overall travel time for a full maximum displacement 
of the density distribution (towards and back). As ar- 
gued before, this time can be obtained from the average 
particle position (z) t by measuring the (time) distance 
between two maxima [see Fig. [3Jb) for an example] . In 
Fig. [3a) we plot the cycle time T as a function of the 



delay time r. Clearly, by increasing the delay time r the 
cycle time T increases as well. To understand this behav- 
ior it is crucial to recall that the time-delayed feedback 
force incorporated here is of Pyragas type, i.e. a control 
target at time t and its value at time t — t is used for 
the feedback signal. In our case the control target is the 
centre of mass position {z) t , which is being shifted as a 
function of time due to the constant force F hlas . Thus, a 
larger delay time r implies that the system travels longer 
distances within the time interval r. On the other hand, 
the specific form of the feedback force F th is constructed 
such that it always counteracts F blas (see Eq. (fTJ) and 
Refs. 0, [HI). Furthermore,^ the absolute value of F ft 
is small for large differences f(f) — f(t — t). As a result, 
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Thus, the initial distribution is given by 



Figure 6. Density oscillation cycle time T (a) as a func- 
tion of the delay time r and (b) as a function of the wash- 
board potential amplitude Uo within the oscillatory regime 
[see Fig.[T]. The parameters are Fo = 7ksT/a, Ko — 7fcsT '/a 
and L — lOna. The rest of the parameters as in the legend. 



the crossover region where the feedback force changes 
from being essentially inactive to compensating the con- 
stant tilting force F blas is accessed more often when the 
delay time is smaller. Thus, smaller delay times r yield 
decreased cycle times T up to the limit where r is too 
small to induce time-periodic oscillations in the density 
any longer [see FigQJb)]. 

Fig. EJb) shows the dependence of T on the washboard 
amplitude Uq. It is seen that the cycle time T decreases 
slightly as a function of the washboard amplitude Uq- It 
is well known that the energy barrier plays a decisive role 
for the particle escape rate in a potential minimum for 
hopping processes that are thermally activated [4l[. In 
our case, the interplay between the washboard potential, 
the constant tilting force and the control force determines 
the rate at which particle are crossing to the next poten- 
tial minimum. By increasing the washboard amplitude 
Uq, the energy barrier for a particle escape is increased 
leading to smaller displacements of the average particle 
position in a valley. As a result the cycle time of the 
oscillations is decreased. 



IV. MICROCHANNEL GEOMETRY 

So far we have considered infinite systems (i.e. systems 
with periodic boundaries). In this section we explore to 
which extent the emergence of an oscillation instability 
also depends on the channel geometry. To this end, we 
now consider a system that consists of a microchannel of 
finite length L with hard walls at the ends. We choose as 
initial distribution the equilibrium density distribution 
for a single particle subject to the washboard potential 
plus a wall potential, j3U waX[ = 10(z/z wa n) 20 , which con- 



< 



^wall 



20(7. 



fines the particle position to values 
Such a smooth wall potential is typically used to model 
situations where the diameter of the particles forming the 
wall is much smaller than that of the fluid particles [42| . 



p(z, t = 0) = po exp 



-PU cos 2 



10 



-wall 



(16) 



Again, we focus on the "balanced case" 7 = Fo — Kq = 
with a fixed control amplitude value of J\o = lk^T/a. 
We recall that in the case of the infinite system (periodic 
boundaries) a washboard amplitude value of Uq = 8fcgT 
is already sufficient to suppress oscillatory states, yield- 
ing a stationary state for times larger than the transient 
response time. We choose this specific value for Uo as a 
starting point for the microchannel study. 

In Fig. [7Ja) we show snapshots of the density distribu- 
tion at three different times to = Ot^, t\ = 64.3tb and 
t 2 = 93.0tb- We note that the initial density distribu- 
tion p(z, t = 0) (shown as a red curve) is symmetric with 
respect to the position z = 0. After a transient response 
time of roughly 20tb we find stable time-periodic density 
oscillations where the distribution oscillates between two 
states characterized by a low and a high average particle 
position, respectively. Therefore, this state fulfills the pe- 
riodicity condition p(z,t + T) = p(z,t) with T = 2t:/uj 
being the cycle time of the oscillations. The snapshots 
at times t\ and t-x = ti + T (shown as a green-dashed 
and a blue-dashed curve in Fig. [TJa), respectively), re- 
veal that, indeed, both density profiles appear to be iden- 
tical, at least to the naked eye. We support this conclu- 
sion by plotting in Fig. [TJb) the average particle position 

(z)t = J Q L zp(z, t)dz as a function of time. Here, the func- 
tion (z)t is averaged over the entire channel length L in 
contrast to the circular ring geometry (periodic bound- 
aries) where the obtained results are periodic with re- 
spect to each valley position. Similar as in Sec. IIII CI 
the function (z) t oscillates as a function of time with 
the same frequency u> as the frequency of the density os- 
cillations, as can be seen from the time stamps t\ and 
ti = t\ + T that we have included as vertical lines. How- 
ever, we stress that the cycle times here are much longer 
than for the circular ring geometry Furthermore, the pe- 
riodic states have the spatial period equal to the largest 
spatial period used in the system, which is the system 
length itself as imposed by the wall potential. As a re- 
sult, the time-periodic solution oscillates back and forth 
between z = ±z wa n. We also note that the finding of 
density oscillations for Uq = 8k gT, 7 = and t = tb 
is in contrast to what we found in Sec. IIII CI where no 
(periodic) instabilities occur for this specific parameter 
set. 

By further increasing the washboard amplitude to the 
value of Uo = 9fc bT, on the other hand, we find that 
the oscillatory behavior at early times is being damped 
resulting in a stationary, non-periodic density for times 
t > 20tb- We conclude that the results of the linear 
stabilty analysis (cf. Fig. [T]) can approximately be used 
as a reference to find states that are linearly unstable 
for the microchannel system. We argue that this is be- 
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Figure 7. (Color online) Results for the microchannel, (a) density distribution p(z,t) as a function of position for selected times 
to/i~B — (red curve), Ii/tb = 64.3 (green-dashed curve) and ti/rB = 93.0 (blue-dashed curve) and (b) shows the first moment 
of p(z,t) as a function of time. The parameters are Uo = 8fcsT, Fo = 7ksT/a, Kq = IksT/a, r = tb and z wa u = ±20u. 
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Figure 8. Density oscillation cycle time T for the microchan- 
nel as (a) function of the delay time r and (b) function of 
the washboard potential amplitude Uo. The parameters are 
F = 7k B T/a, K = 7k B T/a and z wa ii = ±20a. The rest of 
the parameters as in the legend. 



cause the channel is so large (L = 44c) that the system is 
mostly determined by the bulk properties. The instabil- 
ity region for the finite channel seems to be qualitatively 
similar, but increased in size compared to the results for 
the circular ring geometry. 

In Fig. [Ha) we show results for the cycle time T 
as a function of the delay time r for a fixed value of 
A'o = 7fcsT '/a and the choice 7 = Fo — Kq = 0. Holding 
the washboard potential amplitude fixed to the value of 
Uq = Sfc^T, we do not find oscillatory density states be- 
low t = O.Stb- By increasing r, we find a qualitatively 
similar behavior as in the periodic system (see Sec. IIII D[) 
for the cycle time T. Specifically, we find a monotonic 
increase on T as a function of r. On the other hand, 
quantitatively comparing the results for the cycle time 
T to the periodic system reveals that the values increase 



by approximately one order of magnitude. This is clearly 
a consequence of the substantial differences in the oscil- 
lations that we observe for the microchannel: the un- 
derlying nonlinear terms in the dynamical equations [see 
Eqs. (J3ll4j) ] drive the density propagations over the whole 
system size, which is given by the channel length L. Thus, 
it is clear that the cycle time must be significantly longer 
than in the periodic system, where the oscillations oc- 
cur around a valley position with displacements that are 
smaller than na (the spatial period of the washboard po- 
tential) . We note that for the r values that we considered 
(up to t = 10tb), we do not find any upper boundary for 
the linear stability threshold. This behavior seems to be 
similar to the periodic system, where we showed that the 
oscillatory density state cannot be transformed into a sta- 
tionary state by increasing the value of r [cf. Fig. QJb)]. 

In order to investigate the dependence of the cycle time 
T on the washboard amplitude Uq, we hold the delay 
time fixed to the value of r = tb and increase (decrease) 
Uq in steps of AUo = l^sT towards the linear stability 
threshold. We do not find oscillatory density states above 
Uq w 8.63/cbT. Furthermore, we observe a monotonic 
increase of T as a function of Uq, which contrasts with 
the periodic system where we found monotonic decrease. 
Again, this is a consequence of the substantially different 
oscillation mode that we observe for the microchannel. 
As explained above, any oscillatory perturbation applied 
to the system travels over the whole system length. Thus, 
an increased value of Uq now means that the propagation 
of the perturbation is hindered, which results in an in- 
creased cycle time T that corresponds to the travel time 
over a distance of approximately twice the system length 
L (see Scc lnTDll . 

For completeness, we show in Fig. [S] the bifurcation di- 
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Figure 9. Local extrema values of the function (z) t as a func- 
tion of the washboard amplitude Uo- The shaded area marks 
the region where the system exhibits hysteresis. The parame- 
ters are Fo = 7k B T '/a, Ko = 7k B T/a, r = t b , z W aii = ±20a. 



Figure 10. Results for the average particle position (z)t as a 
function of time for different repulsion strengths: eo = 0k B T 
(solid curve), eo = \2k B T (dashed curve) and eo = 15k B T 
(dotted curve). The parameters are Uo = 15k B T, _Fo = 
7k B T/a, K = 7k B T/a, t = t b , L= IOttct and N = 2. 



agram for the microchannel calculated in the same fash- 
ion as in Sec. IIII CI We only find stationary states for 
washboard amplitudes Uq > Uq. Contrary to the infi- 
nite system, however, the average particle position (z) t 
in this stationary state is principally a non-zero constant 
as t — > oo. For simplicity, all stationary states in Fig. [9] 
have been shifted to zero. Upon increase of Uo all density 
oscillations suddenly drop off at Uq = Uq with a jump 
in the (oscillation) amplitude from the value |(z) max | to 
zero; i.e., a subcritical Hopf bifurcation occurs. As the 
parameter Uq is reversed, a stationary solution can be 
found below the Hopf bifurcation point for values rang- 
ing to Uq 1 sa 7 ' .TbksT. Thus, the system system exhibits 
hysteresis within the parameter region U™ < Uq < Uq 
(shown as the shaded area in Fig. [5]). 



V. INFLUENCE OF REPULSIVE PARTICLE 
INTERACTIONS 

In many real colloidal systems, interactions between 
the particles cannot be neglected. Here, we briefly con- 
sider the case of purely repulsive interactions within the 
periodic system. To this end, we utilize the recently devel- 
oped dynamical density functional theory (DDFT). The 
DDFT key equation is given by (43l - |45| 



r -i fy(M) v 

at 



5T[ P {z,t)} 



6p(z,t) 



(17) 



The mobility coefficient in Eq. (jXTJ) is the same as in the 
Fokker-Planck approach [see Eq. Q], i.e. we can set its 
value to Tr^k^T/o- 2 = 1. The chemical potential p(z,t) 
obtained from the Helmholtz free energy functional T 
has three contributions 



(i(z,t) 



6p(z,t) 



Mid CM) + Mint CM) + Mext(M) 

(18) 



The first contribution is the ideal gas term pid = 
fcsTln Ap(z, t) (A denotes the thermal de Broglie wave- 
length), the second contribution p int accounts for particle 
interactions and the third contribution is the external po- 
tential p CKt = U cx t , which in our case includes the contri- 
butions from the tilted washboard potential and the con- 
trol force [see Eq. ©]. In the following, the colloidal in- 
teractions are treated within a mean-field approach, that 
is, Mint (2) = ftz'p(z',t)U rcp (\z - z'\). We employ the 
Gaussian Core Model (GCM) where the interaction po- 
tential is given by £7 rep (|;z — - z'\) — £0 cxp \—{z — z') 2 /cr 2 ] . 
The GCM represents a typical coarse-grained potential 
which describes a wide class of soft macroparticles with 
effective (gyration) radius tr [4(| |43| . We choose positive 
repulsion strengths £0 > 0, such that the interaction is 
purely repulsive. Also, we focus in this section on the cir- 
cular ring geometry (periodic boundaries) and consider 
the case of N = 2. 

We note that, even in the non-interacting case, the 
results from the linear stability analysis (see Sec. IIII Al) 
cannot be used as a basis for comparison. The reason is 
that the calculations in Sec. IIII Al were done with = 1. 
Rather, the washboard amplitude Uq must be increased 
significantly to find a stable stationary state. For exam- 
ple, for the "balanced case" 7 = Fq—Kq = (and £0 = 0) 
we find oscillatory density states for a broad spectrum of 
values for Uq ranging up to Uq ~ Mfc^T (recall that 
Uq — %ksT is sufficient to suppress oscillatory states for 
N = 1 and the rest of the parameters being the same). 
In Fig. [10] we show results for the average particle posi- 
tion (z) t for Uq = lhksT (and A^ = 2). In the case of the 
non-interacting system £0 = (shown as a solid curve in 
Fig. ITUI the function (z) t is indeed being damped as a 
function of time, reflecting a stationary state for times 
larger than the transient response time, which is here of 
the order of t rcs w 2000tb. 

Upon increase of the repulsion strength £q the extrema 
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of the function (z) t increase (see the dashed curve and the 
dotted curve in Fig. [TU1 respectively). This is consistent 
with our earlier finding [27j that repulsive interactions 
support the particles in crossing the barrier, yielding an 
increase of the long-time diffusion coefficient. Moreover, 
for so = 15/cbT we find stable (time-periodic) density 
oscillations with cycle time T = 2.315tb. Thus, repul- 
sive interparticle interactions can be successfully used to 
stabilize the oscillatory density state. We note, however, 
that for the present system increasing the average density 
(via the particle number TV) has a considerably greater 
impact on the linear stability of the system. 



VI. CONCLUDING REMARKS 

In the present work, we have investigated the dynamics 
of colloidal particles subject to a one-dimensional, tilted 
washboard potential under time-delayed feedback control. 
The major goal was to identify conditions under which 
the control can induce oscillatory states. The latter are 
absent in the uncontrolled system. Our investigations are 
based on the (nonlinear) Fokkcr-Planck equation, com- 
bined with a linear stability analysis. 

We have investigated infinite systems (i.e., systems 
with periodic boundaries) and microchannels of finite 
length L (bounded by repulsive walls) . For the first case, 
we have obtained, based on linear stability analysis, a 
full state diagram. This diagram predicts that oscilla- 
tions (of the density field) do indeed occur for finite val- 
ues of the delay time and control strengths comparable 
to the strengths of the conservative forces. Investigating 
the time dependence of the same system via numerical 
solution of the Fokker-Planck equation, we found full con- 
sistency with the results of the linear stability analysis for 
all model parameters considered. In addition, the numer- 
ical solution provides results for the (likewise oscillating) 
moments of the density distribution. In particular, from 
the oscillations of the first moment (i.e., the current) we 
could identify the cycle time. The latter was found to 
monotonically decrease with the delay time. 

In the finite system (which we have investigated via 
the full Fokkcr-Planck equation alone), we also found os- 
cillations. However, the spatial extension of these oscil- 
lations corresponds to the wall-to-wall separation, rather 
than to the width of one valley, as in the infinite-system 
case. Despite these differences, the parameter region 
where the finite system exhibits feedback-induced oscil- 
lations is rather similar to that found in the infinite case. 
We attribute this to the fact that the finite system un- 
der consideration was still so large (L 3> A w b = ttct), 
that boundary effects are not dominant. Finally, we have 
briefly considered the case that the colloidal particles in- 
teract. This was done on the basis of the recently devel- 
oped DDFT (43T - |45| . a generalized continuity equation 
where the particle interactions enter via a free energy 
functional. Using a purely repulsive (GCM) pair poten- 
tial, we have shown that such interactions can stabilize 



oscillatory states. However, to see this effect the strength 
of repulsion must be of the order of the washboard am- 
plitude (and both must be significantly larger than k^T). 
Taken together, we have shown that colloidal particles in 
modulated potentials under time-delayed feedback con- 
trol can display highly non-trivial dynamics, particularly 
oscillations. One way to understand these differences to 
the uncontrolled case is that our feedback control term, 
which relies on the average particle position and thus in- 
volves all particles, introduces effectively time-dependent 
interactions between the particles. 

We note that feedback-induced spatio-temporal effects 
have been recently also been found in other extended 
systems such as optical resonators, where the delayed 
feedback generates spontaneous motion of cavity solu- 
tions (32j, and, more generally, systems describable by 
the Swift-Hohenberg equation [48j| . A particularly inter- 
esting feature of the present colloidal system is that it is, 
in principle, accessible by experiments [l6l - [l8j . Indeed, 
colloidal particles are typically so large that their posi- 
tion (and, consequently, the first moment of the density 
distribution) can be easily monitored by real-space meth- 
ods such as video microscopy. We therefore hope that our 
results will stimulate future experiments. Moreover, from 
the theoretical side one could extend the present analy- 
sis to mixtures consisting of several species with multiple 
time delay constants. This could be a promising route 
for the development of a novel particle sorting effect in 
narrow channels. 
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Appendix A: Numerical continuation of solutions of 
the linearized Fokker-Planck equation 

It is known that in order to continue any given solution 
of the below Sec. IIII Al mentioned boundary value prob- 
lem (BVP) with 7 equations, 7 boundary conditions and 
3 integral conditions, one requires (7 + 3 — 7+1 = 4) con- 
tinuation parameters. We always include the pair {u>, A) 
into the set of continuation parameters. Consequently, 
we are left only with two additional continuation param- 
eters. These can be chosen in the arbitrary fashion out 
of the set (Uq, Kq, Fo,t, N). For example, the solution 
of the BVP can be continued in the plane of (Uq,Ko), 
with all other parameters fixed. The result of such con- 
tinuation is a line in the space (Uq,Kq). Generally, any 
continuation yields a co-dimension one manifold in the 
space of system parameters, which represents the stabil- 
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ity threshold. The onset frequency u is then an output 
parameter, which is given by a certain function of the po- 
sition on the threshold. The accuracy of all continuation 
runs is controlled by checking that the absolute value of 
the eigenvalue A is of the order of 10~ 10 . . . 10~ 12 . 
Generally, it is not easy to find a particular combination 
of the system parameters directly on the stability thresh- 
old, and additionally to correctly guess the corresponding 
value of a;. In order to find a starting point for our con- 
tinuations, we employ the following strategy. First, we 
switch off the integral conditions Eqs. (fTTj) and continue 
an analytically known solution at vanishing driving force 
(F a = 0), that is, 

p.(z)~e- u <'\ C(z) = 0, S(z) = 0, (Al) 



with an arbitrary guess of r and ui in parameter u, until 
the point (S) = is found. Second, we continue this so- 
lution in parameter Kq, until the point on the stability 
threshold is found, i.e. (S) = and (C) = 1. Note that 
the condition (S) = remains uneffected, when continu- 
ing in the parameter Kq. Finally, we switch the integral 
conditions Eqs. (JTTJ back on and proceed with the con- 
tinuation in any of the parameters (Uq, Kq, Fq, t, N), as 
described above. 
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